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Abstract 

CN ■ 

^ . Consider the 3D flow of a viscous Newtonian fluid upon a curved 

2D substrate when the fluid fllm is thin as occurs in many drain- 
ing, coating and biological flows. We derive a comprehensive model 
\^ • of the dynamics of the fllm, the model being expressed in terms of 

^^ ! the fllm thickness rj and the average lateral velocity u. Based upon 

centre manifold theory, we are assured that the model accurately in- 
^ eludes the effects of the curvature of substrate, gravitational body 

^.^. force, fluid inertia and dissipation. The model may be used to resolve 

'^ I wave-like phenomena in the dynamics of viscous fluid flows over ar- 

Q ■ bitrarily curved substrates such as cylinders, tubes and spheres. We 

^ . briefly illustrate its use in simulating drop formation on cylindrical 

flbres, wave transitions, Faraday waves, viscous hydraulic jumps, and 
flow vortices in a compound channel. These models are the most com- 
plete models for thin fllm flow of a Newtonian fluid; many other thin 
fllm models can be obtained by different truncations of the dynamical 
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C^ I equations given herein. 
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1 Introduction 

Mathematical models and numerical simulations for the flow of a thin fllm 
of fluid have important applications in industrial and natural processes p?7| 
!nSl !nHl !nni !Il Ha !2S1 !III!. The dynamics of a thin fluid fllm spreading or 
retracting from the surface of a supporting liquid or solid substrate has long 
been active areas of research because of its impact on many technological 
flelds, for example, coating flows [3^: applications of coating flows range 
from a single decorative layer on packaging, to multiple layer coatings on 
photographic fllm; coated products include adhesive tape, surgical dressings, 
magnetic and optical recording media, lithographic plates paper and fabrics. 
A wide variety of thin fluid fllm models have been reviewed in detail by 
Oron et al. [24J. In this Introduction, we summarise some of the results on 
mathematical models for three dimensional thin fluid fllm flows on a solid 
curved substrate and relate these to the new model derived herein. 

In a three dimensional and very slow flow, a "lubrication" model for the 
evolution of the thickness r/ of a fllm on a general curved substrate was shown 
by Roy et al. [SB] to be 
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1 Introduction 3 

where ( = 1] — ^nrf + |fciA;2'7^ is proportional to the amount of fluid locally 
"above" a small patch of the substrate; k is the mean curvature of the free 
surface of the film due to both substrate and fluid thickness variations {k k. 
K+V^?7); K is the curvature tensor of the substrate; /ci, /c2 and k = ki + k2 are 
the principal curvatures and the mean curvature of the substrate respectively 
(positive curvature is concave); W is a Weber number characterising the 
strength of surface tension; and the differential operator V is defined in a 
coordinate system on the curved substrate. Based upon a systematic analysis 
of the continuity and Navier-Stokes equations for a Newtonian fluid, this 
model accounts for any general curvature of the substrate and that of the 
surface of the film. Deere and Baret ^U] found good agreement between a 
linearised version of lubrication models such as (^ and experiments of flow 
over various shaped depressions in the substrate. 

In many applications the lubrication model (0) of slow flow of a thin 
fluid film has limited usefulness; instead a model expressed in terms of both 
the fluid layer thickness and an overall lateral velocity or momentum flux is 
needed to resolve faster wave-like dynamics in falling films [7, pi 10], wave 
transitions 0, higher Reynolds number flows ^22j Eqn.(19)], in rising film 
flow and a slot coater [13 Eqn.(37)], and in general flows [TUIHII. Roberts j31j 
derived such a model for two dimensional flow, approximately 

dfj d{rju) 



dt dx 
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where TZis a. Reynolds number of the flow, u is the lateral velocity averaged 
over the fluid thickness, and Qqs is the lateral component of gravity. Here 
we greatly extend the model (J2HS1) by deriving in ^the approximate model 
for the flow of a three dimensional thin liquid layer of an incompressible, 
Newtonian fluid over an arbitrary solid, stationary and curved substrate, 
such as the flow about a cylinder shown in Figures Q and El The derived 
accurate model (j^SHSSl) for the film thickness t] and a weighted average lateral 
velocity tt, see (UHl), is to low order and analogous to (I2H3)^ 
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^The approximate equality in the conservation of mass equation Q) becomes exact 
equality when ^ replaces r] on the left-hand side. The higher order analysis leading to (|52|l 
does this automatically. 
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t=0 t=4 



t = 8 t=12 





Figure 1: around a nearly horizontal cylinder of radius R = 2 (not shown) we 
start with a fluid layer of thickness rj = 1 except for a small bump discernible 
off top dead centre of the cylinder. With Reynolds number 7?. = 10 in a 
gravitational field with gravity number Q = 0.5 the fluid bump first slides 
around the cylinder to the bottom by time t = 12 . The other sections of the 
fluid also slide down to the bottom of the cylinder, but not so fast. 

where Qg^ is the component of gravity tangential to the substrate. The con- 
servation of fluid equation (j3]) is the natural generalisation of equation (J2I) to 
three dimensional flow. The momentum equation (jH)) similarly generalises (jHl) 
to three dimensional flow through a nontrivial effect of curvature upon drag. 
The important feature of this model, as in (J2HS1); is the incorporation of the 
dynamics of the inertia of the fluid, represented here by the leading order 
term TZdu/dt , which enables the model to resolve wave-like behaviour. In 
contrast, the lubrication model of thin films (^ only encompasses a much 
more restricted range of dynamics. We base the derivation of model (J3H3) 
upon a centre manifold approach established in ^ The approach is founded 
on viscosity damping all the lateral shear modes of the thin fluid film except 
the shear mode of slowest decay. Then all the physical interactions between 
spatial varying quantities, substrate curvature, surface tension and gravita- 
tional forcing are systematically incorporated into the modelling because the 
centre manifold is made up of the slowly evolving solutions of the Navier- 
Stokes and continuity equations; for example, all these physical effects take 
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t = 24 



t = 36 



t = 48 



t = 60 





Figure 2: around the nearly horizontal cylinder of radius R = 2 (not shown 
but at angle 0.1 radians to the horizontal) the fluid lump now, t = 24, at 
the bottom of the cylinder slowly pulls in fluids from the two ends of the 
cylinder under surface tension of Weber number W = 20 . By times t = 48 
and 60 surface tension forms a large off-centre bead which slowly slides along 
the cylinder, surrounded by a thin layer, rj ^ 0.1 , still covering the cylinder. 



part in the particular simulation shown in Figures G] and 121 For example, the 
7r^/12 coefficient in the models such as ()H|5j] is not 1: the coefficient of these 
terms would be 1 in modelling based upon the heuristic of cross sectional 
averaging; but 7r^/12 = 0.8224 is correct because it must be | of the viscous 
decay rate 7r^/4 in order to match the leading ^ coefficient in the lubrication 
model ((H). In this approach the model is based upon actual solutions of the 
Navier-Stokes equations and so gets all coefficients correct to a controlable 
order of accuracy. 

Before we undertake the construction we introduce in §2] an orthogonal 
curvilinear coordinate system fitted to the substrate. The analysis then starts 
from the incompressible Navier-Stokes equations and boundary conditions 
recorded in ^for this special coordinate system. Thus we derive the model 
for a completely general curving substrate. 

The derived model ()52riH!^ reduces to a model for three dimensional fluid 
flow on flat substrates upon setting the principle curvatures ki and ^2 equal 
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to zero. For example, the low order model (0) becomes simply 



du IT u 



2 „-. ^2 



The higher order accurate version of this model, recorded in §6.11 as ((SHh 
IM|) . extends to three dimensional fluid flows the models for two-dimensional 
fluids on flat substrates that were derived in [SH HI] . In ^6.11 we report on 
the linearised dynamics, t] = 1 + h where both h and u are small. One result 
is that 

7^a,J = -— CU + V^cj, (7) 

where uj = Vx — Uyis a. measure of the mean vorticity normal to the substrate 
in the flow of the film. Thus mean normal vorticity just dissipates due to 
drag and diffusion. However, letting 5 = Ux + Vy , which measures the mean 
divergence of the flow of the film and hence indicates whether the film is 
thinning or thickening, we find 

ht = -5, (8) 

2 2 

nSt = -^5 + ^[ggnV^h + WV^h\+{l + w)V^5, (9) 

where w = 3.0930 . Observe that this divergence diffuses with a larger co- 
efficient, namely 1 + tu , than that of pure molecular diffusion; this effect is 
analogous to the enhanced Trouton viscosity of deforming viscous sheets fIE\ 
pl43, e.g.]. The enhanced viscous dissipation is due to interactions with the 
shear flow similar to those giving rise to enhanced dispersion of a passive 
tracer in pipes [211 e.g.]. From see that the divergence of the film's veloc- 
ity is simply driven by gravity and surface tension acting on variations of the 
film's thickness and is dissipated by substrate drag and the enhanced lateral 
diffusion. Nonlinearities and substrate curvature modify this simple picture 
of the dynamics. 

Circular cylinders are a specific substrate of wide interest. For example, 
Jensen |15j studied the effects of surface tension on a thin liquid layer lin- 
ing the interior of a cylindrical tube and derived a corresponding evolution 
equation. Whereas Atherton & Homsy J2], Kalliadasis & Chang jTHj and 
Kliakhandler et al. ^H] considered coating films flow down vertical fibres and 
gave nonlinear lubrication models. Thus in §6.21 we also record the accurate 
model for flow both inside and outside a cylinder as used in the simulations 
for Figures Hand 121 Axisymmetric flows are often of interest in coating flows. 
Using s as the axial coordinate, to low order a model for axisymmetric film 
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2 The orthogonal curvilinear coordinate system 
flow along a cylinder of radius R is 




(10) 



0.6487^ 



where the upper/lower sign corresponds to flow on the outside/inside surfaces 
of the cylinder. Observe that the curvature of the substrate: modifies the 
expression of conservation of mass; drives a beading effect through rjs/R"^; 
and modifies the drag terms. These are just some special cases of those 
recorded in ^ 

The centre manifold approach we use to derive low-dimensional dynami- 
cal models such as (jlHil) has the advantage that its simple geometric picture 
leads to a complete low-dimensional model ^21] e.g.]. Algebraic techniques 
have been developed for the derivation of a low- dimensional model [33j, the 
correct modelling of initial conditions [011301 EDI and boundary conditions [SI]. 
However, we limit our attention to deriving the basic differential equations of 
the dynamical flow. Other aspects of modelling remain for further study. But 
further, at the end of ^ we investigate high order refinements of the basic 
linearised surface tension driven dynamics of (OQ and determine that the 
model ([52H53|l derived here requires that spatial gradients are significantly 
less than the limit |V?7| < 1.9, for example, see ()57|1 . This quantitative in- 
dication of the extent of the model's spatial resolution is better than that 
for lubrication models such as (^ which require the surface slope to be sig- 
nificantly less than 0.74 instead. Such quantitative estimates of the range of 
applicability are found through the systematic nature of the centre manifold 
approach to modelhng. 



2 The orthogonal curvihnear coordinate sys- 
tem 

In this section, we describe the general differential geometry necessary to 
consider flows in general non-Cartesian geometries. First, we introduce the 
geometry on the substrate, then extend it out into space and establish the 
orthogonal curvilinear coordinate system used to describe the fluid flow. 
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2 The orthogonal curvilinear coordinate system 8 

Let S denote the solid substrate. When S has no unibihcal point, that 
is, there is no point on S at which the two principal curvatures coincide, 
then there are exactly two mutually orthogonal principal directions in the 
tangent plane at every point in S \1'A\ Theorem 10-3]. Let ei and 62 be 
the unit vectors in these principal directions, and 63 the unit normal to the 
substrate to the side of the fluid so that ei, 62 and 63 form a right-handed 
curvilinear orthonormal set of unit vectors. Such a coordinate system is 
called a Darboux frame jlHI. Let xi and X2 be two parameters such that the 
unit tangent vector of a parameter curve X2 =constant is ei, the unit tangent 
vector of a parameter curve a;i=constant is 62, and let y measure the normal 
distance from the substrate. Then on the substrate, points P G 5 , 



with substrate scale factor 



1 dP 
vii dxj 



dP 



TTli 



dXi 



Further, we deduce how the normal varies along the substrate: 

9e3 



dxi 



~ flljtxijCi . 



(12) 
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(14) 



Note that the unit vectors Cj are independent of y. At any point in the fluid, 
written as 

X(xi, X2, y) = P(xi, X2) + yes{xi,X2) , 

the scale factors of the spatial coordinate system are, since positive curvature 
corresponds to a concave coordinate curve. 



hi 



dX 



dxj 
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The spatial derivatives of the curvilinear unit vectors are [3 p598] 
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dxi 



h 



l,V 



hi' 



Bi' + niikies 



-TTliKiCi , 



dy 
dxii 



dy 

hi' A 







hi 



where i' = 3 — i is the complementary index of i and hij denotes dhi/dxj . 

A fundamental quantity for the problem is the free-surface mean curva- 
ture k which is involved in the effects of surface tension through the energy 
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2 The orthogonal curvilinear coordinate system 

stored in the free-surface. As derived by Roy et al. 
curvature of the free-surface 



eqn(37)], the mean 
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where hi = 1 — kirj are the metric coefficients evaluated on the free surface 
(indicated by the tilde), and where 



A = ^hlhl + hlril^ + hhl , 

is proportional to the free-surface area above a patch dxi x dx2 of the sub- 
strate. 

In the analysis we assume the film of fiuid is thin. Conversely, viewing 
this on the scale of its thickness, the viscous fiuid is of large horizontal extent 
on a slowly curving substrate. Thus we treat as small the spatial derivatives 
of the fiuid fiow and the curvatures of the substrate; Deere & Baret |10i pl62] 
report experiments showing that apparently rapid changes in the substrate 
may be treated as slow variations in the mathematical model. Then an 
approximation to the curvature of the free surface is 



K, 



^2 ki 
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+ 0{ti^ + V'^7]) 



1 — kiT] 1 — k2rj 
where the Laplacian is evaluated in the substrate coordinate system as 



(15) 



V^r^ 



17111712 
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dxi \midxij dx2 \m2dx2 



(16) 



For later use, also observe that on the free-surface the unit tangent vectors ti 
and unit normal vector n are 



ti = {hiBi + r]^,e3)/Jhf + i][ 



Ixi -I 



(17) 



fi = (-/i2^xiei - hii]^^e2 + hih2e^)/A . 



'li 



We describe the dynamics of the fiuid using these formula in a coordinate 
system determined by the substrate upon which the film fiows. 
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3 Equations of motion and boundary conditions 10 

3 Equations of motion and boundary condi- 
tions 

Having developed the intrinsic geometry of general three dimensional sur- 
faces, we proceed to record the Navier-Stokes equations and boundary con- 
ditions for a Newtonian fluid in this curvilinear coordinate system. 

Consider the Navier-Stokes equations for an incompressible fluid flow 
moving with velocity field u and with pressure field p. The flow dynam- 
ics are driven by pressure gradients along the substrate which are caused by 
both surface tension forces, coefficient a, the forces varying due to variations 
of the curvature of the free surface of the fluid, and a gravitational accel- 
eration, g, of magnitude g in the direction of the unit vector g. Then the 
continuity and Navier-Stokes equations are 

V-w = 0, (19) 

UU 1 Ui 

^- + uVu = — Vp + -V^u + g. (20) 

at p p 

We adopt a non-dimensionalisation based upon the characteristic thickness 
of the film H, and some characteristic velocity U: for a specific example, 
in a regime where surface tension drives a flow against viscous drag the 
characteristic velocity is U = cr/p and the Weber number W = a /{Up) = 
1 . Reverting to the general case, the reference length is H, the reference 
time H/U, and the reference pressure pU/H. Then the non-dimensional 
fluid equations are 

V-w = 0, (21) 



n 



du ^ 

—-- + u-Vu 
at 



-Vp + V^u + gg , (22) 



where TZ = UHp/p is a Reynolds number characterising the importance of 
the inertial terms compared to viscous dissipation, and Q = gpH'^/{pU) is 
a gravity number analogously measuring the importance of the gravitational 
body force compared to viscous dissipation; the gravity number Q = TZ/J-' for 
Froude number JF = U'^/gH so that when we choose the reference velocity 
to be the shallow water wave speed U = yJgH then JF = 1 and the gravity 
number Q = TZ . 

In the curvilinear coordinate system defined in ^ the non-dimensional 
continuity and Navier-Stokes equations for the velocity field u = uiSi + 
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^262 + f 63 are as follows (adapted from ,3., p599]): 

d d d 

-— (/i2Mi) + 7^— (^1^2) + -7:-{hih2v) = , 
oxi 0x2 oy 
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(24) 
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where u> is the vorticity of the fluid given by the curl 



(jj = 'V X u 
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and where 



^ ui d U2 d d 

These partial differential equations are solved with the following boundary 
conditions. 

1. The fluid does not slip along the stationary substrate, that is 

n = on y = . (25) 



2. The fluid satisfies the free surface kinematic boundary condition 

dr) Ml dr] U2 drj 

7^ = '"~1~7) T7r~ ony = r]. 

ot hi oxi h2 0x2 



(26) 



3. The normal surface stress at the free surface is caused by surface ten- 
sion, in non-dimensional form 



n ■ T ■ n = p + Wk , 



(27) 



where r is the deviatoric stress tensor on free surface, p is the fluid 
pressure at the surface relative to the assumed zero pressure of the 
negligible medium above the fluid, and h is the unit normal to the free 
surface. 
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4 The centre manifold analysis of the dynamics 12 

4. The free surface has zero tangential stress 

t • f • n = , (28) 

where t is any tangent vector to the free surface. We use the two 
hnearly independent tangent vectors in (fTTj) to ensure the boundary 
condition is satisfied for all tangent vectors. 

In this curvilinear coordinate system the components of the nondimensional 
deviatoric stress tensor r are [3, p599] 

- g / 1 ^^» hj^i' ruiki \ 

\hidxi hihi> ' hi ) ' 

1 dv dui rriiki 
hi oxi oy hi 

dv 
rss = 2—. 
oy 

4 The centre manifold analysis of the dynam- 
ics 

We adapt the governing fluid equations ()23H24|1 and the four boundary condi- 
tions (J2SH2HI) to a form suitable for the application of centre manifold theory 
and techniques to generate a low- dimensional dynamical model with firm 
theoretical support. 

Three mathematical tricks place the equations within the centre mani- 
fold framework; these tricks fit the parameter regime of viscous flow varying 
relatively slowly over the substrate. 

1. First we introduce the small parameter e to characterise both the slow 
gradients along the substrate, d/dxi , and the small curvatures of the 
substrate (as curvatures are the partial derivatives of unit normal with 
respect to Xi). e either may be viewed at face value as a mathemati- 
cal artifice or may be viewed as being equivalent to the multiple-scale 
assumption of variations occurring only on a large lateral length scale 
(large compared to the thickness of the fluid). The two viewpoints pro- 
vide exactly the same results. In either case, let the lateral variations 
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scale with the parameter e: 



_d_ 
dxi 



d 

' dx* ' 



ki = ekl , ^2 = ^k. 



2 5 



en 



where * denotes quantities which have been scaled by e. 

2. Second, the presumed small gravitational forcing is treated as a per- 
turbing "nonlinear" effect by introducing the parameter f3 such that 
the gravity number Q = (3"^ . 



3. Third, as introduced by Roberts |^, we also modify the tangential 
stress condition (|^ on free surface, using a parameter 7: at 7 = 
the lateral shear mode of slowest decay actually becomes a marginally 
stable mode; whereas at 7 = 1 the modification vanishes to restore the 
physical stress-free boundary condition ()28|1 . The modification to (jH^ 
is necessary to create the necessary three modes of the centre manifold 
model. Subsequently, evaluating at 7 = 1 removes the modification to 
obtain a model for the physically correct dynamics. 

Now rewrite the governing fluid equations according to the above tricks. For 
convenience, we drop the "*" superscript on all re-scaled variables hereafter. 
Equations (j2nH211) become 



d 
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where the scale factors are h. 
vorticity are 
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■mi{l — ekiy) , and the components of the 
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The boundary conditions ()25H2H|) become 

n = on y = , (32) 

drj ui drj U2 drj 

a7 = '^~^^a ^1~T^ ony = r], (33) 

dt hi oxi /i2 0x2 

7 ~ ~ n mimim2Ui 

ti-T-n = (1-7) on|/ = ?7, (34) 

rjlil 

h- f ■ h = p + Wk, (35) 

where 

h = \Jh1 + e2r/2^ , / = V (e^2??xi)^ + (e/^ir/^J^ + (/ii/i2)2 , 

and the unit tangent vectors ti and unit normal vector fi are 

ti = (hiBi + er]^,^e^)/li , n = {-eh2r]x^ei - ehir]^^e2 + hih2e^)/l . 

The asymptotic expressions for the deviatoric stress f on the free surface are 

Tii = 2e — -— H ^Mj - fcjW + (^(e ) , 

T12 = e — \ Ml U2]+0[t), (36) 

\m2 0x2 rui oxi mim2 mim2 J 

V rrii axi / 



Ti3 = 


dy 


^33 = 


dy 



and the mean curvature of the free surface k, expanded in powers of e, is 
where V^?7 is the same as that in ()16|1 and ^2 = ^i + ^i • 



The tangential stress boundary condition (j34j) has been modified by the 
introduction of the artificial parameter 7. The physically correct boundary 
condition is recovered when 7 = 1. But when 7 = the boundary condi- 
tion (jH^ linearises to 

dui Ui 

-- = — , ony = T], 

dy 7] 
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4 The centre manifold analysis of the dynamics 15 

which gives two neutral horizontal shear modes, Ui oc y . The above equa- 
tions have generalised the physical equations by introducing the extra pa- 
rameters e, 7 and /?. Then by adjoining the trivial equations 

0, (37) 

we obtain a new dynamical system in the variables u, t], p, e, 7 and /3. 
The original system will be recovered by setting e = 1, 7 = 1 and (3 = \/Q . 
However, the two systems are quite different from the view of centre manifold 
theory. The theory now treats all terms that are multiplied by the three 
introduced parameters as nonlinear perturbing effects on the system. So 
the dynamics we describe will be suitable only when there are slow lateral 
variations in Xi of the curvatures of the substrate, oi u, p and 1], small e, and 
a relatively weak gravitational forcing on the system, small Q. In §Hlwe argue 
that evaluating at 7 = 1 is sound, and towards the end of ^Hlwe give evidence 
that lateral variations are slow enough if their logarithmic derivative is less 
than 1.9/77. We now proceed to use centre manifold techniques to develop a 
model of the dynamics. 

A linear picture of the dynamics is fundamental to the application of 
centre manifold techniques to derive a low- dimensional model. The linear 
part of system ()30flHT|) . that is omitting all terms multiplied by a small 
parameter e, 7 or j3, is 

with the boundary conditions (JHIHSSI) linearised to 

u = on y = , 

dri 

V = U on y = rj , 



dt 

duj u. 



on y = 7] , (40) 



dy T] 

dv 

2— p = on y = T] . 

dy 

Note that there are no curvature nor lateral variations in the above linear 
equations, such variations are not included in the leading order approxima- 
tion to the physical system. The linear dynamical system has three types of 
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solutions: first, a motionless film of constant thickness u = v = p = 0,ri = 
constant; second, the family of lateral shear modes Ui ex sin(xj//?7) exp(At) 
where 



A 



X^ 



-JZr]^ 



such that X = tan x ; 



(41) 



and third, the trivial e, 7 and /? being independently constant. Thus, the six 
modes corresponding to zero eigenvalues, the so-called critical modes, are the 
four modes with rj, e, 7 and (3 arbitrarily constant, and the two modes with 
Ui cc y (obtained in the limit x ^ 0). All other modes correspond to negative 
eigenvalues given by ()41|) — they are damped by viscosity. Consequently the 
centre manifold model which we create has six modes: three corresponding to 
critical physical modes; and three corresponding to trivial parameter modes. 

Centre manifold techniques are justifiably applied to infinite dimensional 
dynamical systems whose linearisation has only non-positive eigenvalues and 
such that the nonlinear perturbation terms in the system are smooth and 
bounded, see Gallay [TT]. The perturbation terms in system (jHOHSII), in- 
volving spatial derivatives, are strictly speaking unbounded so the rigorous 
theory does not apply. Nonetheless, by restricting attention to slowly vary- 
ing solutions only, the derivatives remain bounded and the resulting model 
is expected to be accurate [SH]. With this proviso, a low- dimensional model 
of the system may be derived using centre manifold theory. 

Denote the physical fields by v(t) = {ri,ui,U2,v,p) . Centre manifold 
theory guarantees that there exist functions V and G of the critical modes 
where the critical modes evolve in time, that is 



v(t) = V{ri,ui,U2), such that 



d 


V 


dt 


Ui 
. ^2 . 



G(r7,Mi,M2), (42) 



where there is implicit dependence upon the parameters (e, 7, j3) , which are 
treated as small constants, and where Ui are depth-averaged lateral velocities 
defined precisely as 

Ui = - Ui{l - ki,y) dy . (43) 

7] Jo 

This definition ensures that the fluid flux over a point on the substrate is 
simply rju . We proceed to find functions V and G such that v(t) as de- 
scribed by (P^ are actual solutions of the physical equations (|HUH^ satis- 
fying boundary conditions ()32ti35|l . We calculate V and G by an iteration 
using computer algebra |^, see Appendix|Xl In outline, suppose that an ap- 
proximation V and G has been found, and let V' and G' denote corrections 
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we seek to improve V and G . Substituting 



^-^ + ^' I 



V 

Ml 
U2 



G + G' 



into (jHOH^) and its boundary conditions then rearranging, dropping prod- 
ucts of corrections, and using the hnear approximation wherever factors mul- 
tiply corrections (see |HS1 for niore details), we obtain a system of linear 
equations for the corrections. The resulting system of equations is in the 
homological form 

CV' + AG' = R, (44) 

where C is the linear operator on the left hand side of system (jSHHini), ^ 
is a matrix, and R is the residual of the governing equations (jnOHHTj) and 
boundary conditions (jH^HH^ using the reigning approximations V and G . 
The procedure for solving the homological equation P^ is as follows: first, 
choose G' such that R — AG' is in the range of C; second, solve CV' = r.h.s. 
making the solution satisfy the boundary conditions (jHSHSS)) and the defini- 
tions ()43p. Then regard V + V' and G + G' as the new approximation V 
and G respectively. Repeat the iteration until satisfied with the approxima- 
tion. The ultimate purpose is to make the residual -R become zero to the 
required order of error, then the Approximation Theorem [0] in centre man- 
ifold theory assures us that the low-dimensional model has the same order 
of error. The computer algebra program listed in Appendix ^ performs the 
computations. 

5 The high order model of film fiow 

The computer algebra program listed in Appendix Ogives the physical fields 
of slowly varyingthin film fluid flow, and also describes the evolution thereon 
as a set of coupled partial differential equations for the evolution of the film 
thickness rj and the averaged lateral velocities u. 

Computing to low order in the small parameters, it gives the following 
fields in terms of the parameters and a scaled normal coordinate Y = y/rj: 

p = -eWK~e^WV'^ri + eri-^\/ri-u{Y -l){2 + -f/2)-e'^WriK2 
+ eV-u (7(r - 3)/2 - 2{Y + 1)) + ^r/(7„(r - 1) 
+ 0(e3 + # + /5^72); (45) 

u = u{2Y- 7(y^ - Y/2)) + T]^ (e^WVft: + Gg,) [-^lY' 
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240 ' 480 ' 4 2 24 

' V20 ^ 2 ^ 60 ^ 4 ^ 120 ^ 

_ ly3 ^ ^2 li^\ + ^ . ^ (^ A y5 ^ 11 y3 

2 12 / ' VlO^ 60^ 

- ^7>^ -y- \y) + C?(e^ + ^^ + /?^ f ) ; (46) 

t; = eVr/-'u(7(-3r^ + r2)/4 + y2^ 

+ er/V ■ w (7(F^ - F^)/4 - F^) 

+ 0(e3 + M=' + /5^72); (47) 

where Qn is the component of gravity in the direction normal to the substrate. 
The error terms in these expressions involve u = \\u\\ , and then (9(e^ + u'^ + 
/?"*, 7") is used to denote terms s for which either s / {e^ + u'^ + (3^) is bounded 
as (e, -u, /5) — s> , or s/7" is bounded as 7 -^ . The corresponding evolution 
to this order of accuracy is 



dr] 



e V ■ {r]u) + 0{e' + u' + (3\Y), (4^ 



-1 - /3 3 \ 2_ 

- er/ nu ( - - -7 1 - 3r/ u-f 

+ 0{e^ + u^ + P^^''). (49) 

Throughout K is the curvature tensor which in the special coordinate system 
chosen to fit the substrate takes the diagonal form 

Although derived in the special coordinate system, the above and later results 
are all written in a coordinate free form. The differential operators that 
appear are those of the substrate. In the special orthogonal coordinate system 
they involve the substrate scale factors rrii as in the Laplacian ()16|1 . 

To recover the original model, we need to set 7 = 1 so that ()34|) reverts 
to the physically correct boundary condition. In the above asymptotic ex- 
pansions every coefficient is a series in 7, and the ratios of the coefficients of 
7*^"^ to 7" in all such series appear to be greater than about 1.5 for n > 2 
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w/r/2 


9s 


uk/t] 


1 





+0.75000 


-1.50000 


7 


-3.00000 


+0.10000 


+0.60000 


2 

7 


+0.60000 


-0.03286 


-0.10286 


7^ 


-0.06857 


+0.00571 





7^ 





-0.00032 


+0.00321 


7^ 


+0.00128 


-0.00009 


-0.00024 


7^ 


-0.00008 


+0.00002 


-0.00014 


7^ 


-0.00004 


+0.00000 


+0.00003 



Table 1: some higher order terms in the series expansions in 7 of selected 
coefficients in the low-dimensional model (J53|) showing that these expansions 
are effectively summed at 7 = 1 . 



from further calculation. That is, the radii of convergence of the various 
series' in 7 are greater than about 1.5. Table ^ shows the coefficients of 
the 7 series of some terms in a higher order version of the low-dimensional 
model pn|) . Evidently the convergence of at least these series' is very good — 
we may expect five decimal place accuracy from the shown terms and similar 
for the other coefficients. A similar argument on the convergence of such 
series is reported in |311IS2I- Hence we justifiably substitute 7 = 1 into the 
series of every coefficient to obtain the physical model. Here we generally 
calculate every coefficient in the evolution from the terms in the series up 
to and including those of order 7^. We also now set e = 1 to recover the 
unsealed model of the original dynamics. With higher order corrections in 7 
the low-dimensional model ()49|) then becomes the model (0) discussed briefly 
in the Introduction. 



The low order expressions ()45H47|) . when 7 = e = 1 , give approxima- 
tions to the physical state of the fluid flow corresponding to a given t], ui 
and U2 ■ Higher order expressions for the normal structure of the lateral 
velocity u are plotted in Figure El The solid curve shows the fundamental 
structure of the lateral velocity in the normal direction; qualitatively it is 
dominantly parabolic, but it differs in detail. It is indistinguishable from the 
trigonometric | sin(7ry/2) expected from the correct linear problem j32], and 
thus is slightly faster at the free surface than the parabolic proffie with the 
same flux. The dashed curve shows that in order to maintain the flux u the 
flow UiiY) along a trough, fcj/ > , has to be proportionally faster. The dot- 
dashed curve shows that flow curving upwards, around an internal corner, is 
slower at the free surface and conversely faster for flow around an external 
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Figure 3: the normal structure of the lateral velocity field PH|) : solid, com- 
ponent oc u\ dashed, component oc kiiUi\ dot-dashed, component oc kiUi\ and 
dotted, 10 X component oc WVk + Qg^ . 

corner; part of this effect could be attributed to solid body rotation being 
the dissipation free mode for turning a corner. Lastly, the dotted curve, 
exaggerated by a factor of ten, shows the very small adjustment made to 
the profile when the fiow is driven by gravity or lateral pressure gradients — 
observe the velocity at the free-surface will decrease slightly so that when 
lateral forces exactly balance the drag on the substrate the profile will be 
then the familiar parabolic Pouiseille fiow. These show how just some of the 
physical processes affect the details of the physical fields and thus indirectly 
influence the evolution. 

The shear stress on the substrate is of interest: 



'y = 2A67- + 0.1775r] (WVK + gg,) + {Kl- 3.609 K)-u 



+ o{dl + u' + g^/^) 



(51) 



The first term is just the viscous drag on a fiat substrate. The next is 
the enhanced stress transmitted to the substrate when the fiuid is driven 
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by a body force or pressure gradients, equivalently. The third and last term 
accounts for the effects of curvature on the velocity field affecting the velocity 
profile near the bed. 

Higher order model. With computer algebra we readily compute a more 
accurate, more comprehensive higher order model. Atherton & Homsy pOj 
and Lange ^^ have similarly considered high order models of thin film flows 
obtained via computer algebra but only in the lubrication approximation. 
Computing to the next order in spatial gradients e, velocity field u, and 
gravitational forcing (3, the model is written as follows:^ 



dt 

(drag) 

(tension) 

(gravity) 
(advect) 
(viscous) 



—V • {r]u) 



(52) 



TT^ U 



— — + {2K + kI) ■ - 
A if r] 



+ (3.2974 K K - 1.1080 kK + 0.6487 K2I) ■ u 
^2 



+ g 



TX 



--V(fi; + riK2 + V^r/) + 1.0779 rjK -Vk- 0.4891 r]KVK 



12 



[9s + 9n^v) + 0.2554 r^K ■ g, - 0.4891 riKg^ 



+ 



n [1.3464 w ■ Vn + (0.1483 u ■ Vri/r] + 0.1577 V • u) u] 

4.0930_ r^o.4886y . (^0.3461 -> 



nO.8348 



i0.4377 



Vx 



tLv X (.^/^n) 



+ 0.9377-Vr7 x (V x n) - 2.4099^^^ V^ (^^0.8299^ 



(53) 



where the differential operators are those of the substrate coordinate system, 
noting in particular that 3, p599] 



V X u 
V X (63^) 



63- 



1 



ei 



17111712 

1 duj 
m2 8x2 



d{m2U2) d{miUi) 



dxi dx2 



e-2 



TTLi dXi 



^Some of the constants that appear here are tentatively identified: 1.0779 — (tt^ + 
16)/24 , 0.4891 = (tt^ - 4)/12 , 2.4099 = 7rV7 + 1 and perhaps 4.0930 == (Stt^ + 7)/21 . 
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u ■ Vit = Ci 



u ■ Vmi + 



+ 62 




U ■ Vm2 + 

mim2 



Observe that fluid is conserved by (|^. In the above model (j^ for the 
average velocity field we identify the apparent physical source of the terms in 
the various lines by the cryptic words to the left of the lines. Generally the 
viscous drag on the bed, surface tension forces and gravitational forcing show 
some subtle effects of the curvature of the substrate. In faster flows of higher 
Reynolds number, the most usually modelled part of the advection terms, 
the self-advection term u ■ Vu , has the definite coefficient 1.3464 . But note 
that some of the self-advection is also encompassed within the {V ■u)u term. 
This modelling settles, see for example [23 Eqn.(19)] or [SHI EcLn.(23)], the 
correct theoretical value for this and other coefficients. The lateral damping 
via viscosity seems most natural to express in a mixed form involving both 
the general grad-div operator and the curl-curl operator (recall the vector 
identity V^n = V(V ■ w) — V x (V x u)). The involvement of fractional 
powers of the film thickness within the scope of these operators is a convenient 
way to reduce the number of terms within the equation; as yet we have not 
discerned any interesting physical significance to this arrangement. You may 
truncate the above model in a variety of consistent ways depending upon the 
parameter regimes of the application you are considering. 

Lubrication models such as (^ may be derived from (jHSHSS))- Obtain 
simple low order accurate models simply by balancing the drag terms, dom- 
inantly [t^"^ / A)u / rf , against the driving forces expressed by the surface ten- 
sion and gravity terms. This then expresses the average velocity field u as 
a function of the film thickness rj. This is substituted into the conservation 
of fiuid equation ()52|1 to form a lubrication model. Higher order, more so- 
phisticated models are formed by then taking into account the consequent 
time dependence of the weaker previously neglected terms in ()53|) . Any re- 
sultant lubrication model is correct because rational mathematical modelling 
is transitive: a coarser model of a model of some dynamics is the same as 
the coarser model derived directly. 

Slower flow occurs when in a specific application the velocity field is pre- 
dominantly driven by surface tension acting because of curvature gradients, 
whence u = OCVn) . The lateral velocities are significantly damped by vis- 
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cous drag on the substrate. In this case truncate (j3^ to 






TT^ U 



— — + (2K + «/)■- 

4: Iff ^ ' f] 



+ W 



TT 



TT 



VR + 1.0779 r]K -Vk- 0.4891 r/^Vfi: 



Y^(^s + 9n^v) + 0.2554 vK-g,- 0.4891 r^/ty. 



iZ " 



(54) 



That is, you may adjust the dynamical model (J53|l to suit a particular appli- 
cation by choosing an appropriate consistent truncation. 

Order one gradients are encompassed by the model ()53p. 'Long wave' 
models such as (j3^ and (j3^ are based upon the assumption that the lateral 
spatial gradients are small. We here quantify what a 'small gradient' means 
in this context following similar arguments for the Taylor model of shear 
dispersion [2II1 1121 • ^^ modify a simpler version of the computer algebra 
derivation of Appendix ^ to find the centre manifold model of the linear 
dynamics about a stationary constant thickness fluid with surface tension 
but no gravity: 77 = 770 + arj' + C(q;^) and u = au' + 0{a^) where^ 



-f = -770 V ■ tx' + 0(a) 



dr]' 
du' 



(55) 



vl 



,2v72 



.6v76 



8v78 



^0 



-2.47 + 4.09 7/^V^ + .7347/^V^ + .QQllrjlW" + .0223 77=V 



u 



S22 rilV^ + .116r/^V^ + .00168 r/^V^ + .00298 77^V^ 
+ 0(a,V^°). 



7/ 



(56) 



Evidently the linear dynamics are a formal expansion in rj^'V'^ . This ex- 
pansion will converge when the lateral gradients are not too steep. Suppose 
locally a solution has spatial structure approximated by an exponential varia- 
tion, say rj' oc e^^ , then the above expansions become power series expansion 
in 7/qZ/^ . The Domb-Sykes plot (20] of the ratio of successive coefficients in 
FigurelUsuggests that the power series converges for tjqv'^ less than something 
roughly in the range 1/0.15 to 1/0.35. The constant sign of the coefficients 

■^The coefficients in the linear model (|56|) come from evaluating at 7 = 1 an expansion 
with errors 0(7^) • The coefficients used here should be accurate as discussed earlier and 
demonstrated in Tabled 
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Figure 4: Domb-Sykes plot of the formal expansions ()56|) extended to 
0(V^°) , X and +, showing that the radius of convergence, 1/n -^ , may 
be roughly in the range 1/0.15 to 1/0.35 . The © Domb-Sykes plot is for the 
analogous lubrication model ()58|1 showing its radius of convergence is 1/1.83 . 

in ()56|) indicates the convergence limiting singularity occurs for real steep 
gradients. But the strong period 3 oscillations in the Domb-Sykes plot of the 
ratio indicates a complex conjugate pair of singularities occur at an angle of 
±7r/3 to the real axis at nearly the same 'distance'. The generalisation of the 
Domb-Sykes plot to cater for multiple comparable limiting singularities jl21 
indicates that the three singularities are at a distance about 1/0.28 ; that is, 
for any quantity w, the magnitude of the logarithmic derivative of the lateral 
structure 



|V(logu;) 



Vw 



w 



< 



1.9 



(57) 



for the model to converge. For example, apply this limit to the surface thick- 
ness, w = 7], to deduce that the steepness of the fluid variations | V?7| < 1.9 , 
and that accurate approximation is achieved for steepnesses significantly less 
than this rough limit. Hence, steepnesses up to about one should be reason- 
ably well represented by those low order terms appearing in the model (jSHI)- 

For interest, we also investigated the analogous but poorer spatial reso- 
lution of the lubrication model (0) of thin film flow jMl e.g.]. The analogous 
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high order but hnear model is 

^ = W [-.333t]^V^ - .6r7o^V^ - 1.09 t/qV* - 2.00 r^^V^" - 3.67r]^^V^^] t]' 
+ 0(a,V^^). (58) 

Continuing this expansion to errors 0(V^°) see in its Domb-Sykes plot in 
Figure E] that this power series converges only for much less rapid variations 
than the model (|^. For example, the fluid thickness steepness | Vr^l < 0.74 , 
and so should be less than about a third, say, in order for the usual first 
term in the lubrication model rjt = |VVV ■ (r^^V^r/) to be reasonable. Thus 
expect the model (jKH|l developed here to resolve spatial structure roughly 
three times as fine as a lubrication model. 

6 The model on various specific substrates 



The model (J52H53J) is quite sophisticated. It is not obvious how it will ap- 
pear in any particular geometry. Thus in this section we record the model 
for four common substrate shapes: flat, cylindrical, channel and spherical. 
The models are given in terms of elementary derivatives rather than vector 
operators for easier use in specific problems. 



6.1 Flow on a flat substrate resolves a radial hydraulic 
jump 

The simplest example is the flow on a fiat substrate. We: show ID wave 
transitions; simulate Faraday waves; explore divergence and vorticity in the 
linearised dynamics; and lastly demonstrate that modelling the inertia en- 
ables us to resolve hydraulic jumps in a radial flow. 

On a flat substrate use a Cartesian coordinate system [x, y) and let the 
mean lateral velocity u have components u and v respectively (note that in 
this subsection y is a tangential coordinate, not a normal coordinate). The 
substrate has scale factors rrii = m2 = 1 , and curvatures ki = k2 = . The 
model (jn2HSSI) becomes, where Qn is the direction cosine of gravity normal 
to the substrate into the fluid layer and where subscripts on rj denote partial 
derivatives, 

drj_ ^ d{r]u) _ d{r]v) 
dt dx dy 
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TT^ U 



dv 
dy 



n 



du du 

1.5041m— + 1.3464 1;— 

ox oy 



0.1577m 



u 



+ 0.U83- {ur]^ + vr]y) 



+ 



4.0930 



d'^u d'^u 



dx^ dy^ 



3.0930 



dxdy 



+ 4.8333^1^ + ^1^ + 1.9167^1^ + 1.9167^|^ 
Tj ox Tj oy Tj oy rj ox 

^2 



-0.5033^ ^ , nm.i^x n.oo.^- 



+ 0.6094 



7]- 

VyVx 
^2 



, 0.1061^-0.5834- 



- 0.0833^ U 



u 



(60) 



vr^ V 



12 



~^— + -^i^(9y + 9nVy) + W {r]^^y + riyyy)] 



n 



1.3464 m 



dv 
dx 



1.5041m 



dv 
dy 



0.1577m 



du 
dx 



+ 



h 0.1483- (Mr^^ + Mr^j,) 

d'^v d'^v d'^u 

+ 4.0930—^ + 3.0930- 



dx"^ dy^ 

,Vy9v i]^dv 



dxdy 



Vxdu 



h,, du 



+ 4.8333^— + — — + 1.9167—— + 1.9167^ — 



1] dy 



rj dx 

n. 



7] dy 



T] dx 



-0.5033% - ^^ 

1]^ ZT] 



+ 0.1061 






0.5834^ M 



+ (0.6094^ - 0.0833^ 1 M 



Tj^ 



V 



(61) 



Observe that the substrate drag, gravitational and surface tension terms 
are quite straightforward. However, the self-advection terms exhibit subtle 
interactions between the components of the velocity fields due to the specific 
shape of the velocity profiles. Further subtleties occur in the viscous terms 
which not only show explicitly a differential lateral dispersion of momentum, 
but also a complex interaction with variations in the free surface shape. 
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Figure 5: two instants of a 2D fluid falling down a vertical plane substrate 
with Reynolds number 7^ = 20 , gravity and Weber number Q = W = 3 . 
The fluid thickness r^ as a function of distance x shows that white noise at the 
inlet a; = is selectively amplified to solitary pulses that move and merge: 
the close pair of pulses near x ~ 130 in (a) at time t = 145 move and merge 
to the large pulse at a; ~ 165 in (b) at time t = 160 . 



Wave transitions: this model resolves ID wave transitions such as those 
reported by Chang, Demekhin and Saprikin [8 . Their parabolicised Navier- 
Stokes equation (1) corresponds to our non-dimensional Navier-Stokes equa- 
tion ()22j) with Q = W = 3 and our Reynolds number TZ = 156 for Chang 
et al's 6. See the corresponding simulations of our model (jKntiHUjl restricted 
to 2D flow (as in [3^) and forced by white noise at the inlet, although the 
forcing here is a little larger than that of Chang et al. 0: the evolution and 
merger of the solitary pulses are qualitatively the same as that simulated by 
Chang et al. The only noticeable difference is that the retained dissipation 
in our modelling almost entirely removes the surface tension waves in front 
of the solitary pulses. 



Faraday waves Vigorous vertical vibration of a layer of fluid on a flat plate 
leads to a rich repertoire of spatio-temporal dynamics [221123 e.g.], such as 
those shown in Figure El We choose the reference velocity to be the shallow 
water wave speed U = \/gH then the non-dimensional parameters TZ = 
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Figure 6: two instants of unsteady waves on a vertically vibrating flat plate 
simulating the Faraday waves vacillating between ordered patterns (left) and 
irregular patterns (right): for mean fluid thickness 1, no surface tension, 
TZ = Q = 40 , and normal gravity modulated by the factor 1 + 0.55 sin(2.2 1) . 



Q . Achieve the vertical vibration simply by modulating the normal gravity 
in (jnnHS]) by, for example, the factor 1 + 0.55 sin(2.2t) . This frequency is 
roughly twice that of waves with wavelength 5 and see in Figure IHl that these 
waves are generated by a Mathieu-like instability. But involved nonlinear 
interactions lead to complex evolution of the spatial pattern of waves. 



Vorticity and divergence We explore some more of the dynamics de- 
scribed by the model ()59flUT|) . Consider the linearised dynamics of small 
variations on a film of nearly constant thickness: i] = 1 + h{x, y, t) where 
h and the lateral velocity u are small. The linearised versions of (j^nHSJ) are: 



ht 
TZut 

nvt 



-Ux - 
+ (1 + W)U:^3: + M™ + tX^fx- 



y ' 
' 12 



[Qigx + gnhx) + W {h:,^^ + h:,yy)] 



yy 



Jxy 1 



vr^ 7r2 



r^ ^ 12 ^^^^^ ^ ^""^^^ ^ ^ '-^^''^ ^ ^^^^^^ 



+ Vxx + (1 + ■(^)v-uy + ruux. 



yy 
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where w = 3.0930 (the name zu is chosen because its value is coincidentally 
close to vr). 

• First take dy of the second from d^ of the third to deduce equation ((7j) 
governing the mean normal vorticity u = Vx—Uy ■ Observe that linearly 
it is decoupled from the other components of the fluid dynamics: the 
mean normal vorticity simply decays by drag on the substrate and by 
lateral diffusion. 
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Conversely, the first of the hnearised equations together with the di- 
vergence of the second two equations decouple from the mean normal 
vorticity to give (jEHHI) for the film thickness and the mean fiow diver- 
gence 6 = Ux + Vy as discussed briefiy in the Introduction. A little 
analysis shows that in the absence of gravity {Q = 0) this model pre- 
dicts damped waves for lateral wavenumbers 

7l/2 



a > a,. 



-K 



7^W/3 - (1 + tu) 



Numerical solutions of the physical eigenvalue problem agree closely 
with this for TZW > 30 even though the critical wavenumber is as large 
as ttc ~ 0.65 . Recall that the limit (jHTj) on logarithmic derivatives in 
this model implies the wavenumber must be less than 1.9 /rj; here the 
critical Oc ~ 0.65 on a fiuid of depth near 1 is comfortably within the 
limit. Waves cannot be captured by the single mode of a lubrication 
model such as ((H). 

• Lastly, in this linear approximation, lateral components of gravity just 
induce a mean fiow in the direction of the lateral component. 

Substrate curvature and the nonlinear effects of advection and large-scale 
variations in the thickness modify this description of the basic dynamics of 
the fluid film. 



Radial flow with axisymmetry: turn on a tap producing a steady stream 
into a basin with a fiat bottom: you will see the fiow spreads out in a thin 
layer, then at some radial distance it undergoes a hydraulic jump to a thicker 
layer spreading more slowly |^ e.g.]. A model with inertia is essential for 
modelling such a hydraulic jump. Here use polar coordinates {r,6), whence 
the substrate has zero curvature ki = k2 = , but scale factors nii = 1 
and m2 = r . Then describe axisymmetric dynamics by neglecting angular 
fiow and variations while retaining the radial velocity u: 

1 
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(a) thickness and mean velocity 




Figure 7: steady axisymmetric radial flow on a flat substrate: (a) free surface 
thickness t] and mean velocity u versus radius r; (b) streamlines showing a 
recirculation under the hydraulic jump at r ^ 33 . The Reynolds number 
TZ = 15 , gravity number Q = 1 and no surface tension, W = . 



The above equations may be integrated in time to see evolving dynamics 
in a radial flow. Here we focus only upon steady solutions obtained by flxing 
an inlet condition of flow leaving the faucet with prescribed thickness and 
velocity; in Figure [3 the flow has rj = 2.25 and u = 2.62 at radius r = 5, 
and exiting the domain with some prescribed mean velocity at large distance, 
M = 0.16 at r = 50 in Figure [7| Newton iteration then flnds the solutions for 
fluid thickness rj{r) and mean velocity u{r) shown in Figure [Tl^a). See the 
flow spreads out in a thin [r] ^ 1) fast flow before undergoing a hydraulic 
jump at distance r ?^ 33 to a thick (77 ^ 4) slow flow. The streamlines in 
in Figure m^b), obtained from the velocity flelds pHHTTjl . show the presence 
of a recirculation under the jump as also seen in the experiments cited by 
Watanabe, Putkaradze and Bohr [H]. Our model expressed in depth aver- 
aged quantities resolves such nontrivial internal flow structures. 

The steady flow in Figure [7| is near the limit of applicability of the model. 
Although the free surface looks steep in the flgure, the slope is everywhere 



Tony Roberts, February 5, 2008 



6 The model on various specific substrates 



31 



less than 1.08 which, although less than the limit ()57|) identified earlier, 
is about as large as one could reasonably use. For interest, other lateral 
derivatives have the following ranges: !]„ G [— .88, .74] , Ur G [—.15,0] and 
Urr G [—.08, .11] . We also find the hydraulic jump with recirculation in less 
extreme fiows than the one shown in Figure [3 However, we have not yet 
found any fiows with an extra eddy at the surface of the hydraulic jump as 
reported for experiments with larger jumps ^T^ §2.1]. 



6.2 Flow outside a cylinder resolves evolving beads 



Thin film fiow on the outside or the inside of circular cylinders or tubes are 
important in a number of biological and engineering applications. For exam- 
ple, Jensen [TH] studied the effects of surface tension on a thin liquid layer 
lining the interior of a cylindrical and curving tube and derived a correspond- 
ing evolution equation in the lubrication approximation. Our model (JS2HSS1) 
could be used to extend his modelling to flows with inertia. Similarly, Ather- 
ton & Homsy |2], Kalliadasis & Chang jlH] and Kliakhandler et al. ^Hj 
considered coating flow down vertical flbres and similarly derived nonlinear 
lubrication models. Here we record the model (JS2HSS1) as it appears in full 
for flows both inside and outside a circular cylinder. The speciflc model for 
a circular tube which is bent and twisted is left for later work. 

Use a coordinate system with s the axial coordinate and 6 the angular 
coordinate; denote the averaged axial and angular velocities by u with com- 
ponents by u and v, respectively. The substrate has scale factors nii = 1 
and 7712 = R where R is the radius of the cylinder, and curvatures fci = 
and ^2 = T1/-R where the upper/lower sign is for flow outside/inside of the 
cylinder. Then the model on a cylinder is, where here C, = rj ± rf /{2R) , 
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For a nontrivial example, see the beading of fluid on a thin fibre in Figures C] 
and 121 gravity first rather quickly moves a lump of fluid to below the hor- 
izontal cylinder of the fibre; thereafter, surface tension more slowly gathers 
more fluid into the beading fluid which beads both below and above the fibre; 
gravity causes the fluid bead to stabilise off the centre of the fibre (curiously, 
there is more in the bead above the fibre than in the original lump); and 
finally the bead slides along the fibre as it is angled downwards a little. The 
three time scales in this evolution, the fast gravity forced flow, the slower 
surface tension driven flow, and the even longer term sliding, are captured 
in this model. 

Simply obtain axisymmetric flows by setting to zero any derivatives with 
respect to Q, and also setting f7„ = (/e = as non-zero values would break the 
symmetry. The equation for v then just describes the decay of angular flow 
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around the cylinder so also set f = . Thus the axisymmetric model is 
d{riu) 



dC 
dt 



ds 
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4 rj'^ lirj 
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(66) 



recall that the upper/lower sign is for flow outside/inside of the cylinder. 
As in lubrication models i36, p254], see that surface tension on the cylinder 
acts through the term Wrjs/R^ in (jMj) or (jU^ rather like a radially outwards 
body force such as the term QgnVs in ^^■ 



6.3 Flow about a small channel grows vortices 

Consider the flow on a substrate with a small channel aligned downhill. We 
compare this viscous flow with the high Reynolds number experiments of 
Bousmar Hj who modelled turbulent flow over flood plains and channels 
in a flume with water of variable depth but of the order of 5 cm deep. 

First create the coordinate system. Bousmar's channel and flood plain 
had constant shape along the stream, the depth only varied across the stream. 
Thus here let s = xi be the along stream coordinate, r = X2 be the horizontal 
distance across the stream on the substrate,^ and y measure distance normal 
to the substrate. The curved substrate is located a distance d{r) > from 
the sr-plane in a normal direction, see the example d{r) in the middle curve 
of Figure IHl^a). Using, for this subsection, j, k and i as the vertical and two 
horizontal unit vectors, the unit vectors, scale factors and curvature of the 
substrate coordinate system are^ 



k 



1 



e,. 



l + d' 



[i - d'j] 



1 



1 + d'- 



-id'i + j] 



^There is no good reason for using the variable name r for distance horizontally across 
the stream, only that it is next to the letter s in the alphabet. In this subsection the 
variable r is not used to indicate any sort of radius. 



^Useful relationships are: 



I + d' , d" = —krw}^ , and m'^ = —d'krrn^ . 
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m.. 



ks = 0, 



rrir 



rCr 



d" 



(1 + rf'')3/2 



This expression for the curvature kj. is well known. Note: the normal co- 
ordinate y is not the vertical coordinate, and so a flat fluid surface located 
at, say, the location of the reference sr-plane is represented by the varying 
y = d{l + d' )^/^ . Similarly nontrivial, as the channel slopes down at an 
angle -(9 = 0.1 radians to the horizontal and not sideways tilted, the gravita- 
tional forcing is in the direction 



g = smv Eg -\ cos v e^ 



cos i} Br 



(67) 



This coordinate system suits any flow where the substrate is almost arbitrar- 
ily curved in only one direction, not just flow along a channel.^ 

Second, computer algebra gives the model on this substrate as, where 
here C = V ~ K^ff /2 , 
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^The model I^HHTHjl reduces to that for the flow outside/inside a cylinder, (|62II64|I . 
when r — 9 and the substrate scale factors are set to kr — +1/-R and nir — R (only 
the direction of gravity H67|l is incorrect). This algebraic connection occurs despite the 
cylinder not being strictly encompassed by a depth d{r) below any reference plane. Indeed 
the beading flow on a cylinder shown in Figures ^ and |5] was actually obtained using code 
for the model (I68H7()II of this subsection. 
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(70) 



Lastly, as expected, simulations show that fast flow develops in the deeper 
channel and slow flow on the shallow regions, see Figure IHl^a). However, the 
shear in the mean downstream velocity. Figure IHl^a), is unstable to relatively 
weak horizontal vortices that grow in the shear and travel downstream, see 
them in the mean lateral velocity shown in Figure |Hfb);^ analogous notable 
vortices were observed by Bousmar E] in their turbulent flows. As also 
noted by Bousmar, see that the vortices here similarly extend into the shal- 
lows, albeit weakly. 

The vortices apparent in Figure |HI[b) fill the computational domain and thus may 
be an artifice of the domain size. However, otherwise identical simulations on twice the 
channel length show twice as many vortices, whereas simulations in a domain half as long 
again show a rich modulation among vortices of roughly the shown length. We deduce 
that the displayed vortices are not solely an artifice of the computational domain. 



Tony Roberts, February 5, 2008 



6 The 




(a) 



lOOv 
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0.0 

-1.1 




(b) 

Figure 8: (a) an example channel is three times as deep in the middle as the 
surrounding shallows, blue; the substrate has curvature kr, green; and the 
fluid flow is nearly ten times as fast in the channel as in the shallows, red; 
(b) this base flow is unstable to superimposed travelling vortices on the shear 
near the sides of the channel as shown by the cross-channel velocity. Here 
7?. = ^ = 80 on a substrate sloping down at angle ■(9 = 0.1 radians, and with 
no surface tension, W = . 
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The simulation reported here has a change in depth of the substrate 
sufficiently big so that the nonlinear nature of the derived model is certainly 
essential: Deere & Baret ^01 pl55] comment that nonlinear theories are 
needed in viscous flow if the change in substrate profile is bigger than half 
the shallow fluid depth; here the factor is about three, that is, six times the 
linear limit identified by Deere & Baret. 



6.4 Flow on the outside of a sphere 

For flow on the outside of a sphere we use a coordinate system with 6 the co- 
latitude coordinate, (p the azimuthal (longitude) coordinate, and co-latitude 
and azimuthal velocity components u and v respectively. The substrate has 
scale factors rrii = R and m2 = R sin 6 where R is the radius of the sphere, 
and curvatures ki = k2 = —1/R. Note: on a sphere every point is an 
umbilical point; nonetheless, the earlier analysis is valid in this conventional 
spherical coordinate system. Then the model on a sphere is, where here 
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These models look horribly complicated but recall that depending upon the 
application, simpler truncations are often appropriate; two such examples 
are (jSj) and (jH^)- These models have the assurance of centre manifold theory 
that all physical effects are included to the controlable specified accuracy. 



7 Conclusion 

We systematically analysed the Navier-Stokes equations for the flow of a thin 
layer of a Newtonian fluid over an arbitrarily curved substrate. The resulting 
general model (j^SHSHI) resolves the dynamical effects and interactions of in- 
ertia, surface tension, and a gravitational body force as well as the substrate 
curvature. We presented evidence towards the end of ^that this model ap- 
plies to flows where the lateral gradients of the fluid thickness are somewhat 
less than 2, see the more precise limit ()57|) . and (in Q where the time scales 
of the flow are reasonably longer than the decay of the second lateral shear 
mode, that is, longer than 0.045 t^^/j/. The centre manifold paradigm for 
dynamical modelling is based upon actual solutions of the governing Navier- 
Stokes equations, parametrised in terms of cross-layer averages. Further the 
paradigm implicitly arranges the interaction terms between various physical 
processes to support flexible truncation of the model as appropriate for dif- 
ferent parameter regimes; thus the relatively complex model (JH2HSS1) may be 
justifiably simplified as needed by your application. 

To illustrate a range of applications we briefly reported some simulations 
of: wave transitions on a sloping substrate, Faraday waves on a vibrating flat 
plate, and a viscous hydraulic jump in radial flow, see ^6.11 the formation 
and sliding of beads on a cylindrical fibre with surface tension and gravity, 
see ^(i.2( and the generation of vortices in the shear flow between a chan- 
nel and surrounding shallows, see §fJ.H[ These simulations demonstrate the 
resolution of the complex interactions between the varied physical processes 
encompassed by the model. 

Acknowledgment: we thank the Australian Research Council for a grant 
to help support this work. 
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A Computer algebra derives the model 

The REDUCE* program which performs the derivation of centre manifold 
model is listed below. A little explanation is usefully given first. Observe 
that program variables are typeset in teletype font. 

• The physical coordinate system within the program is {xi,X2,y) = 
(x, z, y) with scale factors {hi, /i2, h^) = (hi, h2, 1) , velocity field (mi, M2, v) 
(u, w, v) and pressure p = p ■ Scale factors of the substrate are mi = ml 
and 7712 = ni2 , whereas those evaluated on the free surface are hi = hhl 
and /i2 = hh2 . 

Substrates with specific geometries, as discussed in ^ may be derived 
simply by coding information about their curvatures ki and scale fac- 
tors mj as shown for the three cases discussed in ^Hl 

• Expressions are written in terms of a stretched coordinate system X = 
xi, Z = X2, Y = y/rj, T = t so that the free surface of the fluid film is 
simply Y = 1. In the program we use {X, Z,Y,T) = (xs, zs,ys, ts). 
Then spatio-temporal derivatives transform by the chain rule 



d d Y7]x d 
dxi dX 7] dY' 


d 
dx2 


d 
dZ 


Yvz d 
7] dY 


d d Y7]T d 
dt dT 7] dY' 


and 


d 
dy 


1 d 

7]dY' 


as coded. 









• The amplitudes of the model are {rj,ui,U2) = (h, uu, ww) with h(m,n) 
denoting the spatial derivative qx^qz" ^^^ similarly for uu(m,n) and 
ww(m,n). The evolution of these quantities is given by 

d7] dui du2 

a^ = S^' ^ = S"' ""^ ^ = S"- 

The correctness of the results of this program depend only upon the correct 
coding of the physical equations. The algebraic machinations are repeated 
until the residual of the fluid differential equations and boundary conditions 
are zero to the requisite order. 



®At the time of writing, information about reduce was available from Anthony 
C. Hearn, RAND, Santa Monica, CA 90407-2138, USA, mailto : reduceOrand . org. 
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1 Comment Constructs slowly-varying lubrication model of thin 

2 film 3D fluid flows on arbitrarily curved 2D substrates. 

3 Allows for large changes in film thickness. Incorporates 

4 surface tension, gravity, inertia, slip. The coordinate 

5 system is (x,z,y), velocity (u,w,v), scale factors 

6 (hl,h2,l). 
7 

8 Last updated 16 Sept 2003; 

9 7o improve printing appearance 

10 on div; off allfac; on revpri; 

11 factor d, h,uu, WW, re, gravity, betag,gr, we, kl ,k2,sk; 

12 7, use operator h(m,n) to denote df (h,x,m,z,n) 

13 operator h; 

14 depend h,xs,zs,ts; 

15 let { df (h(~m,~n) ,xs) => h(m+l,n) 

16 ,df (h(~m,~n) ,zs) => h(m,n+l) 

17 ,df (h(~m,~n) ,xs,zs) => h(m+l,n+l) 

18 ,df (h(~m, ~n) ,ts) => df (gh,xs ,m,zs ,n) }; 

19 7, principal curvatures of the substrate 

20 operator kl ; 

21 depend kl,xs,zs; 

22 let { df (kl(~m,~n) ,xs) => kl(m+l,n) 

23 ,df (kl(~m,~n) ,zs) => kl(m,n+l) }; 

24 operator k2; 

25 depend k2,xs,zs; 

26 let { df (k2(~m,~n) ,xs) => k2(m+l,n) 

27 ,df (k2(~m,~n) ,zs) => k2(m,n+l) }; 

28 operator uu; operator ww; 

29 depend uu,xs,zs,ts; 

30 depend ww,xs,zs,ts; 

31 let { df (uu(~m, ~n) ,xs) => uu(m+l,n) 

32 ,df (ww(~m, ~n) ,xs) => ww(m+l,n) 

33 ,df (uu(~m, ~n) ,zs) => uu(m,n+l) 

34 ,df (ww(~m, ~n) ,zs) => ww(m,n+l) 

35 ,df (uu(~m, ~n) ,ts) => df (gu,xs ,m,zs ,n) 

36 ,df (ww(~m, ~n) ,ts) => df (gw,xs ,m,zs ,n) }; 

37 7. use stretched coords: ys=y/h(x,z,t) , xs=x, zs=z, ts=t 

38 depend xs,x,y,z,t; 

39 depend zs,x,y,z,t; 

40 depend ys,x,y,z,t; 

41 depend ts,x,y,z,t; 

42 let { df(~a,x) => df (a,xs)-ys*h(l ,0)/h(0,0) *df (a,ys) 

43 , df(~a,z) => df (a,zs)-ys*h(0,l)/h(0,0)*df (a,ys) 

44 , df(~a,t) => df (a,ts)-ys*gh/h(0,0)*df (a,ys) 

45 , df(~a,y) => df (a,ys)/h(0,0) 

46 }; 

47 7. some abbreviations with appropriate scalings 

48 kx:=d*kl(0,0)$ kz :=d*k2(0,0)$ 

49 eta:=h(0,0)$ etax:=d*h(l ,0)$ etaz:=d*h(0, 1)$ 
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50 7, scale factors of coord system: substrate, fluid & surface 

51 depend ml,xs,zs; 

52 depend m2,xs,zs; 

53 hl:=ml*(l-kx*eta*ys)$ 

54 h2:=m2*(l-kz*eta*ys)$ 

55 hhl:=sub(ys=l,hl)$ 

56 hh2:=sub(ys=l,h2)$ 

57 7o gravitational body force 

58 depend gl,xs,zs; 

59 depend g2,xs,zs; 

60 depend gy,xs,zs; 

61 let { df(gl,xs) => -rh2*df (hl,z)*g2+ml*kl(0,0)*gy 

62 , df(gl,zs) => rhl*df (h2,x)*g2 

63 , df(g2,zs) => -rhl*df (h2,x)*gl+m2*k2(0,0)*gy 

64 , df(g2,xs) => rh2*df (hl,z)*gl 

65 , df(gy,xs) => -ml*kl(0,0)*gl 

66 , df(gy,zs) => -m2*k2(0,0)*g2 

67 }; 

68 7o computes mean across fluid layer 

69 operator mean; linear mean; 

70 let{ mean(ys~~n,ys)=> l/(n+l) 

71 ,mean(ys,ys)=>l/2 

72 ,mean(l,ys)=>l }; 

73 7. solves -df (p,ys)=rhs s.t. sub(ys=l ,p)=0 

74 operator psolv; linear psolv; 

75 let {psolv(ys"~n,ys) => (l-ys~ (n+1) )/(n+l) 

76 ,psolv(ys,ys) => (l-ys~2)/2 

77 ,psolv(l,ys) => (1-ys) }; 

78 7o solves df (u,ys,2)=rhs s.t. sub(ys=0,u)=0 & mean(u)=0 

79 7o and df (v,ys,2)=rhs s.t. sub(ys=0,u)=0 & mean(v)=0 

80 operator usolv; linear usolv; 

81 let {usolv(ys"~n,ys) => (ys~ (n+2)-2*ys/(n+3) )/(n+2)/(n+l) 

82 ,usolv(ys,ys) => (ys~3-2*ys/4)/3/2 

83 ,usolv(l,ys) => (ys~2-2*ys/3)/2 >; 
84 

85 7. may derive model in special geometries if activated 

86 if then 7o flat substrate in Cartesian 

87 let { kl(~p,~q)=>0, ml=>l 

88 , k2(~p,~q)=>0, m2=>l } 

89 else if then 7o cylindrical substrate of radius rad 

90 7o xl axial distance & x2 is angle around 

91 7. sk=+l is flow outside cyl, sk=-l is flow inside cyl 

92 let { kl(~p,~q)=>0 , ml=>l 

93 , k2(0,0)=>-sk/rad , sk~2=>l 

94 , k2(~p,~q)=>0 when p+q>0 , m2=>rad } 

95 else if 1 then begin 7. substrate depth(zs) below reference 

96 7. s=xl longitudinal distance & r=x2 is cross stream location 

97 7. dd=df (depth, x2), m2=\sqrt{l+dd"2}, k2=-dd'/m2"3 

98 depend dd,zs; 
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99 let { kl(~p,~q)=>0, ml=>l 

100 , k2(~p,~q)=>0 when p>0 

101 , df (m2,xs)=>0, df (m2,zs)=>-dd*m2'2*k2(0,0) 

102 , df (dd,zs)=>-m2~3*k2(0,0) 

103 , dd~2=>m2~2-l } end 

104 else if then 7, flat substrate in polar coords 

105 7, xl radial distance & x2 is angle around 

106 let { kl(~p,~q)=>0 , ml=>l 

107 , k2(~p,~q)=>0 , m2=>xs } 

108 else if then 7. outside spherical coords xl=colat x2=long 

109 let { kl(0,0)=>-l/rad 

110 , kl(~p,~q)=>0 when p+q>0 , ml=>rad 

111 , k2(0,0)=>-l/rad 

112 , k2(~p,~q)=>0 when p+q>0 , m2=>rad*sin(xs) } ; 
113 

114 7o centre subspace and null evolution 

115 u:=d*2*ys*uu(0,0) ; w:=d*2*ys*ww(0,0) ; v:=0; p:=0; 

116 gh:=0; gu:=0; gw:=0; 

117 7. and linear approx to various reciprocals 

118 7. rhl, rh2 are the approx of 1/hl, l/h2 

119 rhl:=l/ml$ rh2:=l/m2$ rhhl:=rhl$ rhh2:=rh2$ 
120 

121 7. Use d to count number of uu, ww and derivatives of x & z 

122 7o Throw away this order or higher in uu, ww, d/dx & d/dz 

123 let {d"4=>0, gam'8=>0}; 

124 gravity :=d~2*gr; 7o gravity is the coefficient/order of gravity terms 

125 7. Iterate until converges 

126 it:=0$ 

127 repeat begin 

128 write it:=it+l; 

129 7o approximate mean curvature of the surface 

130 curv:=sub(ys=l,kx+kz +d*rhl*rh2*(df (h2*rhl*etax,x) 

131 +df (hl*rh2*etaz,z)) +(kx"2+kz'2) *eta) ; 

132 7o free-surface deviatoric stress (Batchelor, p600) 

133 txx : =sub (ys=l , 2*d* (rhhl*df (u , x) +rhhl*rhh2*df (hi , z) *w 

134 +rhhl*v*df (hl,y))); 

135 txz:=sub(ys=l,d*(rhh2*df (u,z)+rhhl*df (w,x) 

136 -rhhl*rhh2*(df (hl,z)*u+df (h2,x)*w))) ; 

137 txy:=sub(ys=l,hhl*df (rhl*u,y)+rhhl*d*df (v,x)) ; 

138 tzz : =sub (ys=l , 2*d* (rhh2*df (w , z) +rhhl*rhh2*df (h2 , x) *u 

139 +rhh2*v*df (h2,y))); 

140 tzy:=sub(ys=l,hh2*df (rh2*w,y)+rhh2*d*df (v,z)) ; 

141 tyy:=sub(ys=l,2*df (v,y)) ; 

142 7o omega=curl(q) => del~2q=-curl (omega) (Batchelor p599) 

143 ol:=rh2*(d*df (v,z)-df (h2*w,y)) ; 

144 o2:=rhl*(df (hl*u,y)-d*df (v,x)) ; 

145 oy:=rhl*rh2*(d*df (h2*w,x)-d*df (hl*u,z)) ; 

146 7o vertical momentum & normal surface stress 

147 begin scalar veq,ns; 
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148 veq:=re*( df (v,t)+d*u*rhl*df (v,x)+v*df (v,y) 

149 +d*w*rh2*df (v,z)+d*ml*kl(0,0)*rhl*u'2 

150 +d*m2*k2(0,0)*rh2*w'2 ) 

151 +df (p , y) -gravity*gy+rhl*rh2*d* (df (h2*o2 , x) -df (hl*ol , z) ) ; 

152 ns : = (hh2"2*etax"2*txx+2*hhl*hh2*etax*etaz*txz 

153 -2*hhl*hh2~2*etax*txy -2*hhl~2*hh2*etaz*tzy 

154 +hhl~2*etaz~2*tzz+hhl'2*hh2~2*tyy) 

155 -sub(ys=l ,p) * ( (hh2*etax) ~2+(hhl*etaz) ~2+(hhl*hh2) '2) 

156 -we*curv*((]ili2*etax)'2+(hhl*etaz)~2+(hhl*hh2)'2) ; 

157 ok:=if (veq=0)aiid(ns=0) then 1 else 0; 

158 pd:=eta*psolv(veq,ys)+ns*(rhhl*rhh2)~2; 

159 p:=p+pd; 

160 end; 

161 y. xl -momentum & tangential stress on FS 

162 begin scalar ueq,tsl ,gud,ubed; 

163 ueq:=re*( df (u,t)+d*u*rhl*df (u,x)+v*df (u,y) 

164 +d*w*rh2*df (u , z) +d*w*rhl*rh2* (u*df (hi , z) 

165 -w*df (h2,x))-d*ml*kl(0,0)*rhl*v*u ) 

166 +rhl*d*df (p,x)-gravity*gl+rh2*(d*df (oy,z)-df (h2*o2,y)) ; 

167 tsl:=((hhl"2*hh2-hh2*etax'2)*txy 

168 +hhl*hh2*etax*(tyy-txx)-hhl~2*etaz*txz 

169 -hhl*etax*etaz*tzy )*rhhl'2*rhh2 

170 -(l-gam)*sub(ys=l,u)/eta; 

171 ubed : =sub (ys=0 , -u) ; 

172 ok:=if ok and(ueq=0)and(tsl=0)and(ubed=0) then 1 else 0; 

173 gud: =(-mean(ueq*ys ,ys)-tsl/eta)*3/2; 

174 gu: =gu+gud/re/d; 

175 u: =u+usolv(ueq+2*ys*gud,ys)*eta~2+ubed; 

176 end; 

177 7o x2-momentum & tangential stress on FS 

178 begin scalar weq,ts2,gwd,wbed; 

179 weq:=re*( df (w,t)+d*u*rhl*df (w,x)+v*df (w,y) 

180 +d*w*rh2*df (w,z)+d*u*rhl*rh2*(w*df (h2,x) 

181 -u*df (hl,z))-d*m2*k2(0,0)*rh2*v*w ) 

182 +rh2*d*df (p,z)-gravity*g2+rhl*(df (hl*ol,y)-d*df (oy,x)) ; 

183 ts2 : = ( (hhl*hh2"2-hhl*etaz'2) *tzy 

184 +hhl*hh2*etaz* (tyy-tzz) -hh2*etax*etaz*txy 

185 -hh2~2*etax*txz)*rhhl*rhh2'2 

186 -(l-gam)*sub(ys=l,w)/eta; 

187 wbed : =sub (ys=0 , -w) ; 

188 ok:=if ok and(weq=0)and(ts2=0)and(wbed=0) then 1 else 0; 

189 gwd: =(-mean(weq*ys ,ys)-ts2/eta)*3/2; 

190 gw: =gw+gwd/re/d; 

191 w: =w+usolv(weq+2*ys*gwd,ys)*eta~2+wbed; 

192 end; 

193 y. continuity & bed 

194 begin scalar ceq; 

195 ceq:=-( d*df (h2*u,x)+d*df (hl*w,z)+df (hl*h2*v,y) )/ml/m2; 

196 ok:=if ok ajid(ceq=0) then 1 else 0; 
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197 v: =v+eta*int (ceq,ys) ; 

198 end; 

199 7o uu and ww represent weighted mean velocity 

200 begin scalar ueq,weq; 

201 ueq : =uu (0 , 0) *d-niean(u*h2 , ys) /m2 ; 

202 weq : =ww (0 , 0) *d-mean (w*hl , ys) /ml ; 

203 ok:=if ok and(weq=0)and(ueq=0) then 1 else 0; 

204 u:=u+2*ys*ueq; 

205 w: =w+2*ys*weq; 

206 end; 

207 7o refine approx of hi and h2 

208 hleq:=rhl*hl-l$ rhl : =rhl-hleq/ml ; rhhl :=sub(ys=l ,rhl) ; 

209 h2eq:=rh2*h2-l$ rh2:=rh2-h2eq/m2; rhh2: =sub(ys=l ,rh2) ; 

210 write ok:=if ok and(hleq=0)and(h2eq=0) then 1 else 0; 

211 7. kinematic BC on FS 

212 gh:=sub(ys=l,v-u*rhl*etax-w*rh2*etaz) ; 

213 showtime; 

214 end until ok or(it>19); 
215 

216 if then begin 7o get shear stress on the substrate 

217 txy : =sub (ys=0 , hhl*df (rhl*u , y) +rhhl*d*df (v , x) ) $ 

218 tzy : =sub (ys=0 , hh2*df (rh2*w , y) +rhh2*d*df (v , z) ) $ 

219 gam:=l; 

220 on rounded; 

221 print_precision 4; 

222 coeffn(txy,d,l) 

223 coeffn(tzy,d,l) 

224 coeffn(txy,d,2) 

225 coeffn(tzy,d,2) 

226 coeff(gu,d); 

227 coeff(gw,d); 

228 end; 
229 
230 end; 
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